EXO1 : Enseigner la statistique avec R

AFROMPA’R Bouaké 2025

Auteur·rice

Claude Grasland

Date de publication

2025-03-28

On se propoe de montrer l’intérêt de R pour l’enseignement de la statistique aux étudiants de licence de géographie… Cet exercice correspond typiquement au travail qu’on pourrait donner à des étudiants de licence 2 ou 3 ayant suivi des cours de statistique univariée et bivariée et ayant des connaissances de base en cartographie thématique.

On propose de réaliser une analyse de la distribution de l’Indicateur de Développement humain des régions de Côte d’Ivoire entre 2002 et 2022 à l’aide des données accessibles sur le site Global Data Lab . Le programme est rédigé de telle sorte qu’on puisse facilement refaire l’exercice en prenant d’autres dates (ex. Evolution entre 2012 et 2022) ou d’autres pays d’Afrique de l’Ouest.

Code
library(sf)
library(dplyr)

(1) PREPARATION DES DONNEES STATISTIQUES

  • Consigne : Après avoir chargé les fichiers idh2002.csv et idh2022.csv, construisez un tableau de donnée comportant le code et le nom des régions ainsi que leur population et leur indicateur de développement humain aux deux dates.
Code
# Importe les données
tab1<-read.table("data/exo_stats/idh1992.csv", 
                dec=".",
                sep = ";",
                header=TRUE)
tab2<-read.table("data/exo_stats/idh2022.csv", 
                dec=".", 
                sep=";",
                header=TRUE)

# Sélectionne les variables
tab1 <- tab1 %>% select(code=code, pays=pays, nom=nom, pop1992 = pop, idh1992 = idh)
tab2 <- tab2 %>% select(code=code,  pop2022 = pop, idh2022 = idh)

# Fusionne les deux tableaux
don <- left_join(tab1,tab2)


# Affiche
library(knitr)
kable(don, 
      main = "Tableau de départ")
code pays nom pop1992 idh1992 pop2022 idh2022
CIVr101 CIV Abidjan 2728 0.564 5984 0.569
CIVr102 CIV Yamoussoukro 112 0.623 245 0.619
CIVr103 CIV Bas Sassandra 529 0.406 1160 0.520
CIVr104 CIV Comoe 749 0.378 1642 0.582
CIVr105 CIV Denguele 381 0.332 836 0.440
CIVr106 CIV Goh-Djiboua 686 0.424 1504 0.551
CIVr107 CIV Lacs 1098 0.389 2409 0.542
CIVr108 CIV Lagunes 1023 0.456 2245 0.566
CIVr109 CIV Montagnes 1365 0.380 2994 0.492
CIVr110 CIV Sassandra-Marahoue 1614 0.407 3541 0.512
CIVr111 CIV Savanes 628 0.409 1376 0.471
CIVr112 CIV Vallee du Bandama 520 0.443 1140 0.554
CIVr113 CIV Woroba 574 0.199 1259 0.423
CIVr114 CIV Zanzan 834 0.320 1829 0.491

(2) STATISTIQUE UNIVARIEE

Paramètres principaux

  • Consigne : Etudiez l’évolution de l’IDH entre 1992 et 2022 en vous servant de paramètres principaux (valeurs centrales, paramètres de dispersion). Puis établissez deux histogrammes permettant de visualiser l’évolution.
Code
# sélectionne les variables
sel <- don[,c("idh1992","idh2022")]

# Tableau standard
quant<-apply(sel,2,quantile)
moy<-apply(sel,2,mean)
ect<-apply(sel,2,sd)
cv<-100*ect/moy
tab<-rbind(quant,moy,ect,cv)
row.names(tab) <-c("Minimum","Q1","Médiane","Q3","Maximum","Moyenne","Ecart-type", "C.V. (%)")

kable(tab, caption="Paramètres principaux", digits =2,
      col.names = c("Situation en 1992","Situation en 2022"),
      )
Paramètres principaux
Situation en 1992 Situation en 2022
Minimum 0.20 0.42
Q1 0.38 0.49
Médiane 0.41 0.53
Q3 0.44 0.56
Maximum 0.62 0.62
Moyenne 0.41 0.52
Ecart-type 0.10 0.06
C.V. (%) 24.77 10.65

Histogrammes

  • Consigne : Etablissez deux histogrammes permettant de visualiser la forme de la distribution de l’IDH en 1992 et en 2022
Code
par(mfrow=c(2,1))

mintot <-min(c(sel$idh1992, sel$idh2022))
maxtot <-max(c(sel$idh1992, sel$idh2022))

# Histogramme
hist(sel$idh1992,
     breaks=quantile(sel$idh1992),
     xlim=c(mintot,maxtot),
     col="gray80",
     main= "Situation en 1992",
     xlab = "IDH régional",
     ylab = "Fréquence moyenne")
rug(sel$idh1992, col="black", lwd=2)
lines(density(sel$idh1992), lty=3,lwd=2)

hist(sel$idh2022,
     breaks=quantile(sel$idh2022),
     xlim=c(mintot,maxtot),
     col="gray80",
     main= "Situation en 2022",
     xlab = "IDH régional",
     ylab = "Fréquence moyenne")
rug(sel$idh2022, col="black", lwd=2)
lines(density(sel$idh2022), lty=3,lwd=2)

(3) PREPARATION DES DONNEES GEOMETRIQUES

  • Consigne : Après avoir chargé le fonds de carte affichez le contour des régions et le code des unités.
Code
# Importe le fonds de carte
map <- st_read("data/exo_stats/mapreg.shp")
Reading layer `mapreg' from data source 
  `/Users/claudegrasland1/worldregio/afromapr/data/exo_stats/mapreg.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 14 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -8.599302 ymin: 4.361807 xmax: -2.494888 ymax: 10.73664
Geodetic CRS:  WGS 84
Code
# Projette
map <- st_transform(map, 2042)


# Affichage du fonds de carte
par(mar=c(0,0,3,0))
plot(map$geometry, 
     col="gray90",
     main = "Code des unités spatiales de la zone d'étude")

# Ajout du code des unités spatiales
coo<-st_coordinates(st_centroid(map))
text(coo, map$code, cex=0.5,col="black",)

(4) CARTOGRAPHIE THEMATIQUE

Cartes de stock

  • Consigne : Réalisez deux cartes de stock décrivant le nombre d’habitants en 2002 et 2022 Vous utiliserez la même échelle de taille pour rendre les deux cartes comparables.
Code
library(mapsf)
map<-map[,c("code","geometry")]
map_don <- merge(map, don, by="code")

maxequ<-max(don$pop012,don$pop2022)

par(mfrow=c(1,2))
mf_map(map_don$geometry, col="white")
mf_map(map_don, type="prop", var="pop1992",
       val_max = maxequ, inches=0.1, col="gray20", 
       leg_title = "Population (milliers)",)
mf_layout(title="1992",frame = T, credits = "Source : Global Data Lab, 2025")

mf_map(map_don$geometry, col="white")
mf_map(map_don, type="prop", var="pop2022",
       val_max = maxequ, inches=0.1, col="gray20",
       leg_title = "Population (milliers)")
mf_layout(title="2022",frame = T, credits = "Source : Gloabl Data Lab 2025")

Cartes de ratio (choroplèthes)

  • Consigne : Réalisez deux cartes de taux décrivant le niveau de l’IDH en 2002 et 2022. Pour les rendre comparables vous utiliserez dans chaque carte une partition en quartiles (4 classes d’effectifs égaux)
Code
library(mapsf)
map_don <- merge(map, don, by="code")
maxequ<-max(don$hdi2002,don$hdi2022)

par(mfrow=c(1,2))
mf_map(map_don, type="choro", var="idh1992",
       breaks = "quantile",nbreaks = 4, pal ="Grays",
       leg_title = "IDH",leg_val_rnd = 2)
mf_layout(title="1992",frame = T, credits = "Source : Global Data Lab, 2025")

mf_map(map_don, type="choro", var="idh2022",
       breaks = "quantile",nbreaks = 4, pal ="Grays",
       leg_title = "IDH",leg_val_rnd = 2)
mf_layout(title="2022",frame = T, credits = "Source : Global Data Lab, 2025")

(5) STATISTIQUES BIVARIEES

Nuage de points

  • Consigne : Tracez un nuage de point montrant l’évolution de l’indicateur de développement humain entre les deux dates.
Code
# prépration de l'analyse
code<-don$code
nom<-don$nom
X<-don$idh1992
Y<-don$idh2022
tab<-data.frame(code,nom,X,Y)

# Diagramme
plot(tab$X,tab$Y, 
     pch=20,
     cex=0.8,
     col="red",
     main = "Evolution de l'IDH",
     xlab="IDH 1992",
     ylab ="IDH 2022")
text(tab$X,tab$Y,tab$nom, 
     pos=2,
     cex=0.5,
     col="blue")

Analyse de la corrélation

  • Consigne : calculez les coefficients de corrélation de Pearson et Spearman et testez leur sgnificativité.
Code
cor.test(X,Y, method="pearson")

    Pearson's product-moment correlation

data:  X and Y
t = 4.9223, df = 12, p-value = 0.0003524
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5073463 0.9403486
sample estimates:
      cor 
0.8177876 
Code
cor.test(X,Y, method="spearman")

    Spearman's rank correlation rho

data:  X and Y
S = 134, p-value = 0.006405
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.7054945 

Droite de régression

  • Consigne : calculez l’equation de la droite de régression et tracez- là sur le graphique.
Code
modreg <- lm(Y~X)
summary(modreg)

Call:
lm(formula = Y ~ X)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.052586 -0.016570 -0.001574  0.019274  0.072368 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.33948    0.03848   8.823 1.36e-06 ***
X            0.45013    0.09145   4.922 0.000352 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.03342 on 12 degrees of freedom
Multiple R-squared:  0.6688,    Adjusted R-squared:  0.6412 
F-statistic: 24.23 on 1 and 12 DF,  p-value: 0.0003524
Code
# Diagramme
plot(tab$X,tab$Y, 
     pch=20,
     cex=0.8,
     col="red",
     main = "Evolution de l'IDH",
     xlab="IDH 1992",
     ylab ="IDH 2022")
text(tab$X,tab$Y,tab$nom, 
     pos=2,
     cex=0.5,
     col="blue")

abline(modreg,col="black",lwd=1)

Analyse des résidus

  • Consigne : Calculez les valeurs théoriques prévus par le modèle de régression et les résidus. Affichez le tableau correspondant après l’avoir trié par ordre de résidus croissants.
Code
tab$Y_est <- modreg$fitted.values
tab$Y_res <- modreg$residuals
tab<-tab[order(tab$Y_res),]
kable(tab, digits=3)
code nom X Y Y_est Y_res
11 CIVr111 Savanes 0.409 0.471 0.524 -0.053
5 CIVr105 Denguele 0.332 0.440 0.489 -0.049
1 CIVr101 Abidjan 0.564 0.569 0.593 -0.024
9 CIVr109 Montagnes 0.380 0.492 0.511 -0.019
10 CIVr110 Sassandra-Marahoue 0.407 0.512 0.523 -0.011
13 CIVr113 Woroba 0.199 0.423 0.429 -0.006
3 CIVr103 Bas Sassandra 0.406 0.520 0.522 -0.002
2 CIVr102 Yamoussoukro 0.623 0.619 0.620 -0.001
14 CIVr114 Zanzan 0.320 0.491 0.484 0.007
12 CIVr112 Vallee du Bandama 0.443 0.554 0.539 0.015
6 CIVr106 Goh-Djiboua 0.424 0.551 0.530 0.021
8 CIVr108 Lagunes 0.456 0.566 0.545 0.021
7 CIVr107 Lacs 0.389 0.542 0.515 0.027
4 CIVr104 Comoe 0.378 0.582 0.510 0.072

Cartographie des résidus

  • Consigne : Cartographiez les résidus après les avoir standardisés.
Code
library(mapsf)

# Standardisation des résidus
tab$Y_res_std<-tab$Y_res/sd(tab$Y_res)

# Jointure avec la carte
map<-map[,c("code","geometry")]
map_reg <- merge(map, tab, by="code")

# Choix de la palette et des classes
library(RColorBrewer)
mypal<-brewer.pal(n = 6, name = "RdYlBu")
mybreaks = c(-10, -2,-1,0,1,2,10)

mf_map(map_reg, type="choro", var="Y_res_std",
       pal = mypal, breaks=mybreaks,
       leg_title = "Résidus standardisés",leg_val_rnd = 1)
mf_layout(title="Ecarts à la tendance 1992-2022",frame = T, credits = "Source : Global Data Lab, 2025")

(6) ANALYSE SPATIALE

On se propose d’analyser l’effet de la distance à la capitale économique (Abidjan) ou politique (Yamoussoukro) sur l’évolution de lIDH entre 1992 et 2022

Distance aux deux capitales.

  • Consigne : ajoutez au tableau de données une colonne correspondant à la distance en km à Abidjan ou Yamoussoukro et faites en une cartographie en prenant comme bornes de classes 0, 100, 200, 300, 400, 500, 1000 km.
Code
map_don$capeco<-0
map_don$capeco[1]<-1
capeco<-map_don[map_don$capeco==1,]
map_don$distcapeco<-1+as.numeric(st_distance(st_centroid(map_don),st_centroid(capeco)))/1000

map_don$cappol<-0
map_don$cappol[2]<-1
cappol<-map_don[map_don$cappol==1,]
map_don$distcappol<-1+as.numeric(st_distance(st_centroid(map_don),st_centroid(cappol)))/1000



# Choix de la palette et des classes
mypal<-brewer.pal(n = 7, name = "Greys")
mybreaks = c(0,100,200,300,400,500, 1000)

par(mfrow=c(1,2))
mf_map(map_don, type="choro", var="distcapeco",
       pal = mypal, breaks=mybreaks,
       leg_title = "en km",leg_val_rnd = 0)
mf_layout(title="Distance à la capitale économique",frame = T, credits = "Source : EE CIST 2025")

mf_map(map_don, type="choro", var="distcappol",
       pal = mypal, breaks=mybreaks,
       leg_title = "en km",leg_val_rnd = 0)
mf_layout(title="Distance à la capitale politique",frame = T, credits = "Source : EE CIST 2025")

Relation entre développement et distance à la capitale

  • Consigne : Etudiez les corrélations entre l’IDH en 1992 (Y1) ou 2022 (Y2) et la distance à la capitale économique (X1) ou politique (X2). Que pouvez-vous en conclure ?
Code
X1<-map_don$distcapeco
X2<-map_don$distcappol
Y1<-map_don$idh1992
Y2<-map_don$idh2022

par(mfrow=c(2,2))
plot(X1,Y1,main="Distance à Abidjan 1992", sub = round(cor(X1,Y1),3), pch=20)
abline(lm(Y1~X1),col="red")
plot(X2,Y1, main = "Distance à Yamoussoukro 1992",sub = round(cor(X2,Y1),3),pch=20)
abline(lm(Y1~X2),col="red")
plot(X1,Y2,main="Distance à Abidjan 2022", sub = round(cor(X1,Y2),3), pch=20)
abline(lm(Y2~X1),col="red")
plot(X2,Y2, main = "Distance à Yamoussoukro 2022",sub = round(cor(X2,Y2),3),pch=20)
abline(lm(Y2~X2),col="red")